NASA  CoalPMlor  RcfNiit  IMNS 
ICASE  Report  No.  94-29 


AD-A280  998 


lllllllll 


9  ICASE 

A  THREE  DIMENSIONAL  MULTIGRID 
REYNOLDS-AVERAGED  NAVIER-STOKES 
SOLVER  FOR  UNSTRUCTURED  MESHES 


D  J.  Mavriplis 


L 


DTIC 

ELECTE 

Tt 


Contract  NASI -19480 
May  1994 


QUAUTT  QISFBCIBD  8 


Institote  for  Computer  Applications  in  Science  and  Engine«ing 
NASA  Langley  Research  Center 

Hampton,  VA  23681-0001  ^  ^'0 


94-20407 


Operated  by  Universities  Space  Research  Association 

5  0  62 


j 


A  THREE  DIMENSIONAL  MULTIGRID  REYNOLDS- 
AVERAGED  NAVIER-STOKES  SOLVER  FOR  UNSTRUCTURED 

MESHES 


D.  J.  Mavriplis 

Institute  for  Computer  ^iplications  in  Science  and  Engineoing 
NASA  Langley  Research  Center 
Hampton,  VA 


ABSTRACT 

A  three-dimensional  unstnictured  mesh  Reynolds  averaged  Navier-Stokes  solver  is  described. 
Tuibulence  is  simulated  using  a  single  field-equation  model.  Computational  oveiheads  are 
minimized  through  die  use  of  a  single  edge-based  data-structure,  an  efficient  multigrid  solution 
technique,  and  the  use  of  multi-tasking  on  shared  memory  multi-processors.  The  accuracy  and 
efficiency  of  the  code  are  evaluated  by  computing  two-dimensional  flows  in  three-dimensions 
and  comparing  with  results  from  a  previously  validated  two-dimensional  code  which  employs 
the  same  solution  algorithm.  The  feasibility  of  computing  fliree-dimensional  tuihulent  flows  on 
grids  of  several  million  points  in  less  than  two  hours  of  wall  clock  time  is  demonstrated. 


This  research  was  supported  under  the  National  Aeronautics  and  Space  Administration  under  NASA  Con¬ 
tract  No.  NASl-19480  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications  in  Sci¬ 
ence  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681. 


1.  INTRODUCTION 

In  recent  years,  the  use  of  unstnictuied  meshes  has  become  mme  widespread  for  compu¬ 
tational  fluid  dynamics  problems.  The  advantages  of  unstractured  meshes  lie  in  their  ability  to 
deal  with  aibitrarily  complex  geometries,  while  providing  a  nahiral  setting  for  the  use  of  adap- 
tive  mesh  enrichment  techniques.  However,  most  of  the  successes  of  unstractured  meshes 
have  been  in  solving  inviscid  flows,  particulariy  in  three  dimensions.  Recently,  two- 
dimensional  turbulent  flow  solutions  using  unstractured  meshes  have  been  demonstrated 
[1, 2,3,4].  However,  few  attmpts  at  solving  similar  three  dimensional  flows  are  known.  While 
many  publications  have  appeared  in  the  literature  concemirtg  three-dimoisional  unstractured 
grid  Navier-Stokes  computations,  few  if  any  of  these  have  demonstrated  engineering  quality 
solutions  for  practical  aerodynamic  flows.  The  inclusion  of  viscous  terms  into  an  existing 
unstructured  Euler  solver  is,  in  fact,  a  simple  proposition.  The  real  challenge  is  to  devise  a 
complete  solution  strategy  capable  of  resolving  complex  flows  with  good  accuracy  at  accqrt- 
able  cost  on  highly  stretched  grids. 

The  problems  associated  with  unstracmred  mesh  computations  of  turbulent  viscous  flows 
ate  flueefold.  Firstly,  a  suitable  mesh  with  highly  stretched  elements  in  the  boundary  layer  and 
wake  regions  must  be  generated.  Secondly,  a  turbulence  model  capable  of  operating  efficiently 
on  unstractured  meshes  must  be  incorporated.  Finally,  the  memory  and  CPU  overheads  associ¬ 
ated  with  the  solutitm  tedmique  must  be  low  enough  to  allow  for  the  use  of  very  fine  meshes 
which  are  required  for  meaningfol  results.  Much  work  has  been  performed  in  the  two- 
dimensional  setting  in  order  to  alleviate  these  problems.  However,  the  extension  of  these  ideas 
to  three  dimensions  has  often  been  hindered  by  the  overwhelming  computatimial  overheads 
incurred  by  most  mettKMlologies. 

This  p^r  describes  the  development  of  an  efficient  three  dimensional  Navier-Stokes 
solver.  Most  of  the  techniques  employed  have  been  developed  and  draionstrated  previously  in 
the  two-dimensional  setting,  and  this  work  involves  their  extension  to  three  dimensions.  By  a 
careful  choice  of  data-stractures,  the  use  of  a  rapidly  converging  multigrid  algorithm,  and  the 
implementation  of  parallel  processing  techniques,  a  solution  technique  which  incurs  acceptable 
overheads  and  is  capable  of  dealing  with  relatively  fine  grids  is  obtained. 

2.  METHOD  DESCRIPTION 


2.1.  Discretization  and  Data  Structures 

We  seek  steady-state  solutions  to  the  Favre-averaged  Navier-Stokes  equations.  These 
equations  must  be  closed  using  a  suitable  turbulence  model  in  order  to  model  the  Reynolds- 
stress  terms.  Spatial  discretization  of  the  Navier-Stokes  equations  is  performed  using  a  Galer- 
kin  finite-element  approach.  The  conserved  flow  variables  are  stored  at  the  vertices  of  the 
mesh,  and  the  convective  fluxes  are  assumed  to  vary  linearly  over  each  tetrahedral  element. 

For  the  viscous  terms,  velocities  are  assumed  to  vary  linearly  over  the  tetrahedral  elements. 
Velocity  gradients  can  thus  be  constructed  at  element  centers,  which  then  enables  the  discreti¬ 
zation  of  the  second  derivatives  contained  in  the  viscous  terms.  Additional  artificial  dissipation 
terms  are  constructed  as  a  blend  of  Laplacian  and  biharmonic  operators.  The  Laplacian  dissipa¬ 
tion  results  in  locally  first-order  accuracy,  and  is  thus  triggered  only  in  the  vicinity  of  shock - - 

waves,  while  the  third-order  accurate  biharmonic  dissipation  is  employed  throughout  the 
flowfield. 
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The  natural  data-stnicture  which  arises  from  a  finite-element  point  of  view  is  the  element 
data-structure,  in  which  a  list  of  elements  is  stored,  with  pointers  for  each  elemoit  identifying 
the  four  vertices  whidi  constitute  that  element  It  has  previously  been  shown,  in  the  context  of 
inviscid  flow  calculations,  that  the  convective  terms  can  be  assonbled  using  an  edge-based 
data-structure,  which  is  both  more  compact  (in  terms  of  memory  overheads),  and  minimizes  the 
amount  of  gather-scatter  requited  on  vector  and  parallel  computer  architectures.  The  basic 
data-stracture  for  assonbling  the  convective  terms  is  thus  a  list  of  edges.  Fbr  each  edge  we 
store  the  addresses  of  the  two  end-points  of  the  edge,  and  three  coefficients,  which  represent 
the  x,y,and  z  components  of  the  normal  of  the  face  of  foe  dual  mesh  pierced  by  foe  edge,  as 
shown  in  Figure  1,  for  the  two-dimensiorud  case.  In  three  dimensions,  if  one  considers  all 
tetrahedra  which  share  a  givoi  edge,  as  shown  in  Hgure  2,  the  face  normal  associated  with  this 
edge  can  be  comimted  as  the  sum  of  all  foe  tetrahedral  faces  which  touch  only  one  of  foe  two 
end-points  of  the  edge  (i.e.  sum  of  all  the  in  Figure  2). 

When  employing  a  Galerkin  finite-element  discretizalion,  the  viscous  terms  are  tradition¬ 
ally  foought  of  as  a  sequence  of  two  loops,  one  to  construct  gradients  at  triangle  or  tetrahedron 
centers,  and  another  to  form  the  final  residual  contributions.  However,  the  firud  discrete 
viscous  terms  obtained  in  this  maruier  form  a  nearest  neighbor  stencil.  The  viscous  terms  for  a 
vertex  i  depend  only  on  values  at  i  and  at  vertices  k,  such  that  k  is  joined  to  i  by  a  mesh  edge. 
Thus,  an  edge-based  data-structure  may  also  be  employed  to  assemble  the  viscous  terms.  This 
fact  has  previously  been  pointed  out  in  several  references  [2,5,6].  In  [2],  a  complete  derivation 
of  the  edge-based  coefficients  for  a  Hessian  matrix  is  given.  In  three-dimensions,  fois  would 
require  the  storage  of  9  coefficients  per  edge,  sitKe  the  discrete  Hessian  is  written  as 
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where  Vol„  represents  the  volume  of  the  union  of  tetrahedra  which  touch  vertex  i.  However, 
the  local  edge-based  coefficient  matrix  is  synunetric  about  the  diagonal.  This  fortunate  fact  has 
apparoitly  gotK  uruioticed  in  the  literature.  A  proof  of  fois  is  given  in  the  Appendix.  Thus, 
we  need  only  store  6  coefficients  per  edge  for  the  discretization  of  the  viscous  terms.  If  the 
cross  derivative  terms  in  foe  Navier-Stokes  equations  are  neglected,  fois  can  be  reduced  to  3 
coefficients  per  edge.  Also  note  that  for  foe  discretization  of  a  Laplacian,  a  single  coefficient 
per  edge  is  required,  which  is  given  by  the  sum  of  the  diagonai  terms.  Thus,  if  we  adopt  the 
thin-layer  form  of  the  Navier-Stokes  equations,  a  sit^le  edge  coefficient  is  sufficient  to  com¬ 
pute  the  viscous  terms. 


At  this  stage,  the  foil  Navier-Stokes  terms  have  been  included,  and  thus  6  edge- 
coefficients  are  stored.  It  is  useful  to  examine  the  amount  of  storage  this  represents  compared 
to  other  iqrproaches.  Tetrahedral  unstructured  meshes  of  N  points  contain  aN  tetrahedra  and 
(a  +  1)N  edges,  neglecting  boundary  effects.  The  value  of  a  depends  on  the  mesh,  but  usually 
lies  between  5  and  6.  The  traditional  elemoit  data-structure  requires  the  storage  of  foe  four 
comers  of  each  tetrahedron.  Taking  a  ==6,  fiiis  corresponds  to  24N.  Excessive  storage  is 
required  if  a  douUe  loop  is  employed,  with  foe  cell  gradients  stored  as  intermediate  values  (at 
least  9  extra  values  per  cell  or  54N).  A  single  loop  over  the  tetrahedra  can  however  be  used  to 
assemble  the  foil  viscous  terms.  Within  the  loop,  the  gradients  may  be  computed,  divided  by 
the  volume,  multiplied  by  the  appropriate  face  normal,  and  then  accumulated  to  one  of  the  four 
comer  vertices  of  the  cell.  The  last  two  operations  are  repeated  for  each  comer  vertex  of  the 
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cell.  This  approadi  requires  only  ttie  storage  of  the  etement  dtoa-stnicture  {dus  the  vertex 
coordinates  (which  are  not  required  in  the  edge-based  approach),  thus  a  total  of  27N.  The 
edge-based  api»oach  requires  a  total  of  42N  or  ISN  more  than  die  elonent-based  apfwoadi. 
However,  the  elonent  approach  requires  an  (Mder  of  magnitude  more  operations  since  the 
geometric  quantities  such  as  cell  volume  and  face  normals  must  be  recomputed  each  time  in 
the  loop.  If  any  of  these  quantities  are  stored  rather  than  recomputed,  the  storage  requirements 
quickly  exceed  those  of  the  edge-based  approach.  This  is  a  classic  example  of  a  memoiy-cpu- 
dme  trade-off. 

Since  the  viscous  terms  of  the  Navier-Stokes  equations  are  not  stricUy  second  derivatives 
of  a  given  quantity,  but  contain  terms  such  as 

V.  ( nV,tt  )  (2) 


an  exact  implanentadon  of  a  Galerkin  discretization  requires  the  evaluation  of  the  viscosity  at 
the  tetrahedra  centroids  in  the  construction  of  the  edge  coefficients,  (see  ^ipendix).  If  the 
viscosity  is  not  constant,  the  edge  coefficients  caiuiot  be  computed  in  a  preprocessing  phase. 
Thus  some  approximation  must  be  made  in  order  to  oiable  the  use  of  preprocessed  edge 
co^^cients.  The  simplest  approach  is  to  iqiproximate  the  viscosity  for  eadi  edge  contribution 
as  the  average  of  the  two  values  at  each  end  of  the  edge.  Thus  the  viscous  terms  are  computed 
as 


(PUxlx 

dm,), 

(im,). 

(P«,), 

dm,). 

dm,). 

(P«.)x 

dm,). 

dm,). 

1 

Vol„ 


imthbon 

z 


* 

Oyx 

Ou 

ft 


On 

“yy  “y* 

0.1 


(  M.  +  Pt  ) 
2 


iii  -  «*  ]  (3) 


It  is  important  to  realize  that  this  discretization  does  not  corre^nd  exactly  to  the  Galerkin 
discretization.  Only  under  ceitain  conditions  do  the  two  formulations  become  equivalent  (e.g 
linear  viscosity  variations  and  regular  triangulations).  Another  approach  is  to  rewrite  the 
viscous  terms  as 


V  ( pVx  )  =  nV^x  +  VpVx  (4) 

and  evaluate  the  double  gradient  term  at  each  vertex  using  a  nearest  neighbor  stencil.  Gradient 
calculations  of  this  type  can  be  performed  using  the  convective  edge-based  coefficients.  Both 
of  these  approaches  have  been  tested  in  two-dimensions  and  virtually  no  difference  in  the  final 
solution  could  be  seen.  These  modifications  are  necessary  for  the  viscosity  in  the  momentum 
equations,  the  values  of  \i*velocity  in  the  energy  equation,  as  well  as  for  the  eddy  viscosity  in 
the  turbulence  equation.  The  effect  of  the  eddy  viscosity  variation  is  expected  by  far  to  be  the 
most  important  of  all  for  practical  turbulent  flows.  In  the  present  work,  the  formulation  of 
equation  (3)  is  employed  exclusively. 

Finally,  another  approach  which  has  been  suggested  for  computing  the  viscous  terms  [6], 
consists  of  forming  the  gradients  of  velocity  at  each  vertex,  using  the  nearest  neighbor  stencil 
of  each  vertex,  and  forming  the  second  derivatives  of  the  viscous  terms  by  reapplying  the  same 
integration  to  tire  first  derivatives.  This  approach  is  attractive  because  it  requires  little  addi¬ 
tional  storage.  All  operations  can  be  performed  using  the  edge-based  coefficients  required  for 
the  convective  terms.  However,  it  can  easily  be  seen  that,  on  a  one  dimensional  mesh  of  spac¬ 
ing  h,  this  scheme  reduces  to  a  second  difference  on  a  stencil  of  size  2h.  This  will  result  in 
lower  accuracy  and  possible  odd-even  decoupling.  Since  packing  enough  points  into  the 
viscous  layers  is  generally  one  of  the  main  difficulties  associated  witii  viscous  flow 
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oompiiialioiis,  a  adwme  that  openles  on  evety  other  pdnt  is  highly  underirhUe.  Tlds  scheme 
should  theiefiDie  be  rejected. 

2J.  TurbdaMC  Modeifaig 

The  one-eqaati(Mi  tuibulence  model  of  Spalart  and  Alfanaras  [7]  has  been  im[demenled  in 
the  presem  solvw,  since  dtis  t^jpitMcfa  avoids  die  com|dicalions  involved  in  imidementing  alge¬ 
braic  modds  on  unstructured  meshes  [8],  is  leasrauMy  robust  and  inerqiensive,  and  has  been 
shown  to  yield  favondrie  results  in  three  dimensional  aerodyiuunic  flows  [9]. 

The  tufbulenoe  equation  is  also  discretized  using  the  edge-based  data-stnicture.  The  crai- 
vective  terms  are  treated  using  a  first-order  upwind  fcmnulation,  and  the  difiusirai  terms  are 
treated  analogously  to  die  viscous  terms  in  die  flow  equations.  The  particular  discretization 
ensures  positivity  of  the  discrete  solution.  The  turbulence  equation  is  advanced  in  time  using  a 
point-implicit  treatment  [10].  For  a  single  equation,  this  corresponds  to  a  Jacobi  iteration. 
This  sdieme  is  attractive  since  it  ensures  positivity  of  the  turbulence  equation  variaUe  not  only 
at  steady-state,  but  diroughout  the  convergence  process.  The  turbulence  equation  is  effectively 
solved  decoupled  from  the  flow  equations  (using  diffemit  local  time-step  values).  The  conver- 
goice  of  both  the  turbulence  equation  and  the  flow  field  equations  is  accelerated  using  an 
unstructured  multigrid  algorithm.  The  turbulence  and  flow  equations  are  only  coupled  on  the 
finest  grid  of  the  multigrid  sequence,  where  the  eddy  viscosity  values  from  the  turbulence 
model  ate  fed  back  into  the  flow  solution. 

This  particular  field-equation  turbulence  model  requires  information  concerning  the  dis¬ 
tance  of  each  grid  point  to  the  closest  wall  boundary.  The  simplest  way  of  computing  this  dis¬ 
tance  function  is  to  compute  tlK  normal  distance  from  a  given  grid  point  to  all  boundary  faces, 
and  preserve  the  minimum  distance  found.  This  results  in  an  algorithm,  where  N 

represents  the  number  of  grid  points.  For  fine  grids,  this  becomes  excessively  costly. 
Although  more  so(Msticated  search  techniques  exist  for  reducing  this  to  0{NlogN),  a  simple 
solution  is  to  compute  the  distance  function  on  the  second  finest  mesh  of  the  multigrid 
sequence,  and  then  to  interpolate  these  values  to  the  finest  grid.  This  alone  reduces  the  cost  of 
computing  the  distance  fimction  by  a  factor  of  32. 

23.  Solution  Technique 

The  flow  equations  and  turbulence  equations  are  advanced  in  time  to  obtain  file  steady 
state  solutioa  The  flow  equations  are  advanced  using  a  multi-stage  Runge-Kutta  explicit 
scheme,  while  the  turbulence  equation  is  solved  as  described  above.  Local  time-ste[ping  and 
residual  averaging  are  employed  to  accelerate  the  convergence  of  the  flow  equations. 

An  unstructured  multigrid  technique  is  employed  to  further  accelerate  the  convergence  of 
both  the  flow  and  turbulence  equations.  This  technique,  which  has  previously  been  demon¬ 
strated  for  the  Navier-Stokes  equations  in  two  dimensions  [1],  and  for  the  Euler  equations  in 
three  dimensions  [5],  employs  a  set  of  non-nested  coarse  and  fine  meshes.  In  a  preprocessing 
step,  the  indirection  arrays  which  correspond  to  the  restriction  and  prolongation  operators 
(interpolation  of  variables,  residuals,  and  corrections)  between  each  successive  pair  of  grids  ate 
constructed,  using  an  efficient  search  algorithm.  On  domain  boundaries,  which  may  not  coin¬ 
cide  between  the  various  grids  due  to  the  discretizaticm  of  curved  surfaces  which  constitute  the 
boundaries,  the  search  and  interpolation  procedures  ate  carried  out  in  the  parametric  space 
which  defines  the  surface  patches  of  the  boundary  geometry.  Linear  interpolation  is  used  to 
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transfer  flow  variables,  residuals  and  corrections  b^ween  the  various  meshes  of  die  multigrid 
sequence.  This  requires  the  storage  of  four  addresses  and  four  coefficieitts  per  mesh  point 
Once  these  operators  have  been  cmistructed  and  stored,  the  inter-mesh  multigrid  transfers  can 
be  implemented  as  a  simple  gather-scatter  of  array  elemoits  within  each  multigrid  cycle. 

The  unstructured  multigrid  algorithm  incurs  ap(MOximately  a  30%  memory  oveihead, 
mostly  due  to  the  extra  storage  required  for  the  coarse  grids  of  the  multigrid  sequoice,  and 
requires  roughly  90%  more  CPU  time  per  cycle,  but  results  in  conveigence  rales  which  are  an 
order  of  magnitude  hi^r  than  foe  single  grid  solver. 

2.4.  Parallel  Processing 

The  use  of  one  basic  and  simple  data-stiucture,  and  foe  choice  of  an  explicit  schone  aug¬ 
mented  with  multigrid,  enable  a  relatively  simple  and  efficioit  parallel  implementation  of  foe 
solver  (Ml  shared  and  distributed  memory  parallel  architectures.  Previous  experience  with  dis¬ 
tributed  memory  machines  has  been  less  than  satisfactory.  Such  implementations  recjuire  con¬ 
siderable  effort  and  have  not  been  able  to  demonstrate  siqierior  performance  for  this  class  of 
problems  [11,12].  The  shared  memory  vector-parallel  architecture  of  the  CRAY-YMP-C90  has 
proved  to  be  very  effective  for  foe  current  type  of  problems.  FOr  example,  the  16  processor 
CRAY-YMP-C90  was  found  to  provide  more  than  double  foe  performance  of  the  INTEL  Delta 
machine  using  S12  processors  on  a  three-dimensional  unstructured  multigrid  Euler  solver  [11]. 
Furthermore,  foe  shared  memory  architecture  and  the  parallel  compiler  support  available  on 
CRAY  machines  enable  a  relatively  simple  implemeruation  using  standard  FORTRAN  77. 
Parallelization  on  shared  memory  architectures  is  imperative  for  large  cases,  not  only  to  speed 
up  solution  time,  but  also  to  avoid  idling  processors  in  a  time  sharing  environment  when  jobs 
becomes  large  enough  to  fill  foe  main  memory  of  the  machine. 

The  majority  of  foe  work  in  foe  present  solver  involves  loops  over  edges.  These  loops 
must  be  both  vectorized  and  parallelized.  Since  multiple  edges  meet  at  each  mesh  vertex,  foe 
loops  contain  data  dependencies  which  inhibit  both  vectorization  and  parallelization.  In  order 
to  vectorize  these  loops,  foe  list  of  edges  is  split  into  subgroups,  or  colors,  such  that  within 
each  color,  no  vertex  dependencies  exist.  The  overall  l(X)p  is  hence  transformed  into  an  outer 
loop  over  all  colors,  and  an  inner  vectorizable  within  each  color.  A  simple  parallelization  stra¬ 
tegy  is  to  further  divide  the  colored  groups  into  subgroups  that  can  be  computed  in  parallel. 
This  is  automatically  done  at  compile  time  by  foe  autotasking  compiler  provided  foe  appropri¬ 
ate  compiler  directive  is  specified  at  foe  begirming  of  each  loop.  The  subgroups  are  then  distri¬ 
buted  over  all  processors,  taking  advantage  of  foe  complete  vector  and  parallel  power  of  the 
machine. 

For  the  inviscid  version  of  foe  present  solver,  speedups  of  13  to  14  on  16  processors 
have  been  observed  in  a  dedicated  environment,  indicating  a  degree  of  parallelism  of  over  99% 
has  been  achieved  according  to  Amdahl’s  law.  The  viscous  cases  presented  in  the  results  sec¬ 
tion,  however,  were  run  in  a  time  sharing  environment,  and  yielded  speedups  between  10  and 
12  on  16  processors.  Benchmarking  of  foe  full  viscous  solver  in  a  dedicated  environment  is 
plaimed  for  the  near  future. 


3.  RESULTS 
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3.1.  Sin^  Segmart  Wing 

In  order  to  validate  the  three  dimensional  Navier-Shdces  solver  a  two-dinmnsicHUl  tur¬ 
bulent  flow  over  a  wing  geometry  with  no  qianwise  variatkm  has  been  computed,  and  the 
results  compared  with  the  solution  finxn  a  two-dimaisional  unstructured  flow  solver  on  an 
equivaloit  grid.  The  particular  two-dimensional  solver  employed  has  been  previously 
described  and  validated,  and  is  routinely  used  in  production  envinMunents  [1,13].  This  two- 
dimensional  solver  emfdoys  the  equivalent  discretization,  solution  tedmitpie,  and  mrbuloice 
model  as  the  three-dimensional  code  described  in  this  paper.  Therefore  the  two  codes  should 
give  nearly  identical  results  for  purely  two-dimensional  cases. 

The  first  test  case  involves  the  transonic  flow  over  a  wing  of  aspect  ratio  2  with  no 
sweep  or  spanwise  variation.  The  wing  section  (independent  of  span  location)  is  an  RAE  2822 
airfoil.  The  three-dimaisional  grids  employed  for  computing  the  flow  over  this  wing  geometry 
are  displayed  in  Figures  3.  They  are  formed  by  first  constructing  a  two-dimensional  unstruc¬ 
tured  grid  about  an  RAE  2822  airfoil,  using  the  method  described  in  [14].  The  two- 
dimensional  mesh  is  then  stacked  in  the  spanwise  direction,  thus  forming  a  mesh  of  spanwise 
prizms.  This  prizmatic  mesh  is  then  converted  into  a  tetrahedral  mesh  by  dividing  each  prizm 
into  three  tetrahedra  using  a  variant  of  the  prizm  division  algorithms  reported  in  [15,16]  .  The 
resulting  geometry  consists  of  a  wing  with  a  symmetry  plane  at  both  ends  of  the  wing.  There 
is  thus  no  wing  tip  present  and  no  spanwise  variation  whatsoever.  This  can  be  thought  of  as  a 
typical  wing-in-wind-tunnel  two-dimensional  test. 

The  finest  mesh  for  this  case  is  depicted  in  Figure  3a.  The  entire  mesh  contains  1.04 
million  points  and  6  million  tetrahedra.  The  mesh  is  formed  by  33  spanwise  stations  with 
31,571  grid  points  at  each  station.  The  normal  mesh  spacing  at  the  wing  surface  is  KT^ 
chords,  which  results  in  cell  aspect  ratios  of  the  order  of  500:1  in  these  regions.  A  total  of 
five  meshes  were  used  in  the  multigrid  sequence  for  this  case.  Four  of  these  meshes  are  dep¬ 
icted  in  Figures  3a  through  3d.  Each  coarser  mesh  contains  a  factor  of  approximately  8  fewer 
points  than  the  previous  mesh,  and  consecutive  meshes  are  generally  non-nested. 

The  freestream  Mach  number  for  this  case  is  0.73,  the  incidence  is  2.79  degrees,  and  the 
Reynolds  number  is  6.5  million.  The  solutitm  is  depicted  qualitatively  in  Hgure  4,  as  a  plot  of 
density  contours  on  the  surface  of  the  wing  and  synunetry  walls.  The  lack  of  any  spanwise 
variation  of  the  contours  on  tiie  wing  indicate  the  presence  of  purely  two-dimensional  flow. 
The  flow  is  transonic  and  a  normal  shock  is  observed  slightly  aft  of  the  nud  chord  locatioit 
Hgure  5  provides  a  more  quantitative  picture  of  tiiis  solution.  The  computed  surface  pressure 
at  tile  mid-qian  location  is  compared  witii  experimental  data  as  well  as  witii  the  computed 
results  of  the  two-dimaisional  code  on  an  equivalent  station  grid  of  31,571  points.  Both  the 
two  and  three-dimensional  codes  tend  to  underpredict  the  lift  compared  with  the  experimental 
data,  a  fact  which  is  attributed  to  tiie  turbulence  model  employed.  For  example,  the  same 
two-dimensional  code  adiieves  a  lift  value  some  10%  higher  using  the  Baldwin-Lomax  model. 
However,  die  two  and  three-dimensional  flow  solutions  agree  very  well  with  each  other.  The 
three-dimensional  solution  is  slightly  more  diffusive  than  the  two-dimensional  solution,  which 
is  attributed  to  the  presence  of  extra  spanwise  dissipation,  which  is  non-zero  even  in  a  two- 
dimensional  flow,  due  to  the  presence  of  diagonal  edges  in  between  neighboring  spanwise  sta¬ 
tions. 
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The  conveigence  ntte  for  this  case  is  (dotted  in  Hgure  6.  The  residiials  aie  seen  to  be 
reduced  by  almost  4  orders  of  magnitude  over  200  mutdgiid  cycles,  for  an  average  residual 
reductton  of  0.957.  This  is  com{)arable  but  somewhat  slower  than  the  rate  of  0.943  adiieved 
by  the  two-dimensional  multigrid  code,  which  is  also  (dotted  for  comparison  purposes  in  the 
same  figure.  Fiiudly,  the  single  grid  convergence  rate  of  the  two-dimensi(mal  code  is  also 
shown  in  the  figure.  A  resklual  reduction  of  only  2  orders  of  m^nitude  over  1000  cycles  is 
achieved  for  the  single  grid  approach.  Thus  the  multigrid  procedure  converges  over  10  times 
faster  than  the  single  grid  code,  for  the  two-dimoisional  case.  Since  the  ttiree-dimaisional 
multigrid  convergence  rate  is  close  to  that  of  two-dimensional  case,  one  can  conclude  that 
gains  of  similar  magnitude  are  afforded  by  the  multigrid  algorithm  in  three  dimensions.  The 
three-dimensional  single  grid  convergence  is  not  plotted  due  to  the  excessive  com[mter  costs 
required  for  such  a  run  with  obvious  conclusions. 

In  order  to  demonstrate  the  three-dimensional  capability  of  the  present  code,  the  wing  of 
the  previous  case  (aspect  ratio  2)  is  given  a  sweep  of  30  degrees,  and  a  taper  of  0.5.  This 
spanwise  variation  results  in  folly  three-dimensional  flow.  The  flow  over  the  swept  and 
tapered  wing  at  the  same  conditions  is  com{)uted  using  equivalent  grids  to  foose  described 
above  (same  point  densities).  Figure  7  de(Hcts  the  solution  in  terms  of  density  contours  on  the 
wing  and  symmetry  wall  surfaces.  The  conveigence  rate  for  this  case  is  almost  identical  to 
that  displayed  in  Figure  6  for  the  unswept  three-dimensional  wing,  and  is  therefore  not 
displayed. 

For  both  of  these  cases  a  total  of  177  Mwoids  of  memory  was  required.  This  translates 
to  170  words  per  fine-grid  vertex.  This  includes  all  the  arrays  for  the  coarse  grid  variables  of 
the  multigrid  sequence.  Both  of  these  cases  required  75  seconds  of  CPU  time  per  multigrid 
cycle,  or  a  total  time  of  4.2  single  CPU  hours  on  the  CRAY-YMP-C90  machine.  When  run  in 
multitasking  mode  on  all  16  processors  of  the  machine,  the  total  aggregate  CPU  time  rose  to  5 
hours,  but  the  jobs  were  executed  in  28  minutes  of  connect  time,  as  indicated  by  the  batch  job 
log  file,  in  a  time-sharing  environment  in  which  only  60%  of  tiie  machine  was  dedicated  to  this 
particular  job.  In  a  dedicated  environment,  these  jobs  can  be  expected  to  execute  in  just  over 
15  minutes. 

3J.  Three-EleriMnt  High-Lift  Wing 

The  next  test  case  consists  of  flow  over  a  high-lift  wing  configuration  with  a  slat  and  a 
single  slotted  flap.  The  wing  has  an  aspect  ratio  of  2  with  no  sweep  or  spanwise  variatioa 
The  wing  section  (independent  of  span  location)  is  a  Douglas  three-element  airfoil,  which  has 
been  extensively  tested  both  numerically  and  experimentally  [13].  The  three-dimensional  fine 
grid  employed  for  computing  the  flow  over  this  wing  geometry  is  di^layed  in  Figure  8.  This 
grid  is  formed  by  first  constructing  a  two-dimensional  unstructured  grid  about  the  three-element 
airfoil  using  the  method  described  in  [14],  and  then  stacking  the  grids  and  subdividing  the 
resulting  prizmatic  elements,  as  described  previously.  The  resulting  geometry  can  be  thought 
of  as  a  typical  wing-in-wind-tunnel  two-dimensional  test,  with  symmetry  end  walls  at  both 
extremities  of  the  wing. 

The  finest  mesh  for  this  case  contains  1.84  million  points  and  10.6  million  tetrahedra. 
The  mesh  is  formed  by  33  spanwise  stations  with  55,865  grid  points  at  each  station.  The  nor¬ 
mal  mesh  spacing  at  the  wing  surface  is  10^  chords,  which  results  in  cell  aspect  ratios  greater 
than  1000:1  in  these  regions.  A  total  of  five  meshes  were  used  in  the  multigrid  sequence  for 
this  case.  Each  coarser  mesh  contains  a  factor  of  approximately  8  fewer  points  than  the 
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previous  me^,  and  consecutive  meshes  are  generally  non-nested.  This  level  of  mesh  resolu- 
tkm  corresponds  to  the  eq^valertt  minimum  required  resolution  for  adequate  performance  pred¬ 
iction  using  the  two-dimensional  code,  as  determined  a  grid  resolution  study  [13].  To  the 
author’s  kiwwiedge,  diis  r^nesents  the  largest  aerodynamic  unstructured  grid  computation 
attempted  to  date. 

The  fieestream  Mach  number  for  this  case  is  0.2.  the  incidence  is  16.21  degrees,  arKl  the 
Reynolds  number  is  9  million.  The  flow  is  assumed  to  be  fully  turbulent,  thus  no  transition 
points  are  specified.  The  solution  is  depicted  qualitatively  in  Figure  9.  as  a  plot  of  density 
contours  on  tl^  surface  of  the  wing  and  Mach  contours  on  the  surface  of  the  synunetry  wall. 
The  lack  of  any  spanwise  variation  of  the  contours  on  the  wing  indicates  the  presence  of 
purely  two-dimensional  flow.  The  computed  surface  pressure  at  the  mid-span  location  is  com¬ 
pared  with  experimental  data  as  well  as  with  the  computed  results  of  the  two-dimensional  code 
on  an  equivalent  station  grid  of  55,865  points  in  Figure  10.  Excellent  agreement  between  the 
experimental  data  and  the  two-  and  three-dimensional  codes  is  observed. 

In  order  to  demonstrate  the  three-dimensional  cap^liQr  of  the  present  code,  the  seg¬ 
mented  wing  of  the  previous  case  is  given  a  sweep  of  30  degrees,  and  a  taper  of  0.5.  This 
spanwise  variation  results  in  fully  three-dimensional  flow.  The  fireestream  Mach  number  is  0.2 
as  previously,  but  the  incidence  is  lowered  to  12  degrees.  The  flow  over  the  swept  and  tapered 
wing  (aspect  ratio  2)  is  computed  using  equivalent  grids  to  those  described  above  (same  point 
densities).  Figure  11  depicts  the  solution  in  terms  of  density  contours  on  the  wing  surfaces 
and  Mach  contours  on  the  wall  surface.  The  wall  Mach  contours  are  qualitatively  similar  to 
those  displayed  in  Figure  9  for  the  previous  case,  while  spanwise  variation  of  the  wing  surface 
contours  is  noted,  indicating  the  presence  of  fully  three  dimensional  flow. 

In  Figure  12,  the  multigrid  convergence  rates  of  the  two-dimensional  code,  and  the 
three-dimensional  unswept-wing  and  swept-wing  runs  are  compared.  The  two-dimensional  run 
and  the  three-dimensional  unswept  wing  run  converge  at  very  similar  rates,  achieving  a  resi¬ 
dual  reduction  of  4.5  to  5  orders  of  magnitude  over  3(X)  cycles.  The  three-dimensimial  swept 
wing  run  converges  somewhat  slower  than  the  previous  two  runs,  but  still  achieves  a  similar 
level  of  residual  reduction  over  450  cycles.  This  similarity  between  the  two-dimensional  and 
three-dimensional  convergence  rates  is  a  good  indication  that  the  full  benefit  of  the  multigrid 
algorithm  has  been  achieved  for  three-dimensional  flows. 

Both  of  file  above  cases  required  a  total  of  312  Mwords  of  memory,  and  140  to  f75 
seconds  of  CPU  time  per  multigrid  cycle  on  the  CRAY-YMP-C90  machine,  depending  on  the 
amount  of  concurroicy  achieved  during  a  particular  run  in  the  time-shaiing  envirorunoit.  A 
typical  450  cycle  run  on  all  16  processors  required  1.9  hours  of  connect  time,  as  indicated  by 
the  batch  job  log  file,  in  a  time-faring  environment  in  which  65%  of  the  machine  was  dedi¬ 
cated  to  this  particular  job.  In  a  dedicated  environment,  such  a  job  can  be  expected  to  execute 
in  approximately  1.25  hours. 

4.  CONCLUSIONS 

The  aim  of  this  paper  has  been  to  demonstrate  the  feasibility  of  computing  turbulent 
viscous  flows  over  complex  geometries  on  fine  grids  using  the  unstructured  mesh  approach. 
While  it  is  recognized  that  this  capability  is  expensive,  it  is  entirely  feasible  from  a  technical 
standpoint.  This  is  made  possible  by  several  factors: 

a)  the  reduction  of  memory  overheads  through  the  use  of  a  single  efficient  data-structure 

b)  the  implementation  of  a  rapidly  converging  multigrid  algorithm 
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c)  the  me  of  parallel  {Mocessing 

d)  the  availability  of  a  laige  central  memory  madiine  with  multiple  rapid  processors  and  simide 
to  use  parallelization  tools. 

The  code  currently  requires  approximately  170  words  per  6ne  grid  vertex,  which  includes 
all  the  coarse  grid  information.  It  is  estimated  that  this  can  be  further  reduced  by  10  to  1S%. 
The  code  runs  at  320  Mflops  on  a  single  CRAY-YMP-C90  processor  and  requires  72 
microseconds  per  vertex  per  multigrid  cycle,  with  200  to  400  cycles  usually  required  for  con¬ 
vergence.  Predicted  performance  on  a  dedicated  16  processor  machine  is  4  Gflops  and  S.8 
microseconds  per  vertex  per  cycle,  corresponding  to  a  speedup  of  13.  Thus,  on  the  current 
maximum  configuration  of  the  machine  (1  Gword,  16  processors),  a  case  of  6  million  grid 
points  could  be  computed  in  2  to  4  hours.  Of  course  whether  such  a  computation  is  worth  its 
cost  in  an  industrial  environment  is  another  matter. 

One  of  the  main  obstacles  to  employing  this  capability  in  a  production  environment 
remains  the  grid  generation  process.  The  ability  to  reliably  generate  very  fine  highly  stretched 
three-dimensional  tetrahedral  grids  is  still  under  development,  and  recent  progress  has  been 
made  in  this  area  [15,16].  At  this  stage,  the  flow  solver  has  been  demonstrated  by  constmcting 
three-dimensional  unstructured  grids  using  a  two-dimensional  stacking  procedure.  In  the 
future,  solutions  over  more  complex  geometries  such  as  partial  flap  wings  and  complete  aircraft 
configurations  will  be  attempted.  This  will  require  interfacing  with  a  more  general  grid  genera¬ 
tion  procedure. 
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APPENDK  A 


SYMMETRIC  EDGE-BASED  COEITICIENTS  FOR  VISCOUS  TERMS 


1.  2D  Case 

The  edge  coefficients  for  the  viscous  tenns  are  derived  by  considering  the  Galerkin 
finite-element  discretization  of  these  terms.  Assuming  piecewise  linear  variation  of  the  veloci¬ 
ties,  the  Navier-Stokes  equations  are  multiplied  by  a  test  function  of  compact  support  and 
integrated  over  the  entire  domain.  This  derivation,  which  is  detailed  in  [1],  yields  the  follow¬ 
ing  discrete  form  for  the  second  derivatives  at  vertex  i: 


(P«x)> 

3  y 

f-(p«,)  Ay, 

(pa,)  Ax, 

2Vol,,  t 

[-(H«j,)  Ay, 

(P«,)  Ax, 

where  the  summation  is  over  all  triangular  elements  which  have  a  vertex  at  i.  Vol„  represents 
the  volume  or  surface  covered  by  the  union  of  these  triangles,  and  Ax,  ,Ay,  denote  die  x  and  y 
increments  along  the  outer  edge  of  the  triangle  e,  as  depicted  in  Figure  1.  The  first  derivatives 
represent  gradients  on  the  triangles  of  the  mesh.  For  piecewise  linear  functions,  such  gradients 
can  be  evaluated  exactly  : 


where  the  cell  volume  is  given  by 


(2) 


vol^  =  AxitAyij  -  AxijAyn 


(3) 


for  triangle  A  shown  in  Figure  1.  Using  the  properties  of  differences  such  as  Axy  =  -Axji  and 
Axij  +  Axjt  =  Axu ,  the  volume  can  be  rewritten  as 

vol^  =  AXji^Ayn  -  AxuAyjt  (4) 


By  analogy,  the  gradients  and  volume  for  triangle  B  are  given  by 


vols  =  AxiiAyjt  -  AxflAyu 


(5) 


In  evaluating  the  coefficient  for  the  edge  joining  vertices  i  and  j,  it  can  be  seen  that  only  trian¬ 
gles  A  and  B  in  the  summation  will  yield  non-zero  contributions.  Substituting  equations  (2) 
and  (5)  into  (1),  and  setting  (Ax,, Ay,)  =  (Axy*4yy*)  for  element  A,  and  (Ax,  Ay,)  =  (Axij,Ay,j)  for 
element  B,  we  obtain  a  final  expression  for  the  double  derivatives.  For  example,  flie  xx 
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derivative  at  vertex  i  takes  the  form 


VolJ\luA^  =  I 

^  tdgts 


^yut  ^yjt 


^yu  ^yji 
voIb 


(6) 


where  the  summation  is  over  all  incoming  edges  for  vertex  i.  The  viscosity  values  and  \Ib 
represent  cell  average  values.  Since  this  involves  the  values  at  nodes  I  and  k  as  well  as  at  i 
and  j,  this  construction  is  inconsistent  with  a  purely  edge-based  formulation.  For  this  reason, 
these  values  are  ^proximated  by  the  average  of  the  values  at  vertices  i  and  j,  denoted  as  Pov- 
The  four  double  derivatives  can  now  be  written  as: 
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+ 
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The  terms  involving  grid  metrics  can  all  be  precomputed  and  stored  as  edge-based  coefficients. 
Furthermore,  the  metric  expressions  are  symmetric  about  the  i-j  indices.  Therefore  the  same 
coefficients  may  be  used  for  constructing  the  flux  contribution  to  vertex  i  from  vertex  j,  as  the 
flux  contribution  to  vertex  j  fiom  vertex  i. 

The  remaining  property  to  be  proved  is  the  equality  of  the  two  cross-derivative 
coefficients  i.e 


voIa 


^yu  ^JI  ^  Axa  Ayjt  ^  Ax,;  Ay,-, 
voIa 


(8) 


voIb  voIa  voIb 

This  property  is  not  evident  at  first  glance.  Certainly  each  term  on  the  left  hand  side  of  (8)  is 
not  equal  to  the  corresponding  term  on  the  right-hand  side.  In  order  to  prove  the  equality,  we 
form  the  difference  of  the  two  expressions  in  (8),  and  insert  the  expressions  for  the  cell 
volumes  in  the  denominator  taken  from  (4)  and  (5): 

^yik  Axyt  Ay,,  Axj, 


Axy*  Ay,*  -  Axu  Ay;* 
Ax,*  Ayy* 


Ax,,  Ayy,  -  Ax;,  Ay„ 
Ax,;  Ay;, 


(9) 


Ax;*  Ay,*  -  Ax,*  Ay;*  Ax,,  Ay;,  -  Ax;,  Ay,, 
which  can  be  rearranged  as  follows 

Ax;*  Aya  -  Axg  Ay;*  _  Ax,;  Ay;,  -  Ax;,  Ay,, 

Ax;*  Ay^  -  AXtt  Ay;*  Ax,,  Ay;,  -  Ax;,  Ay,, 

Hence,  the  equality  of  the  two  cross-derivative  coefficients. 


=  0 


(10) 
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2.  3D  Case 


An  analogous  constniction  of  the  edge-based  coefficients  in  3D  yields  the  expression 


= 


(11) 


where  the  correspond  to  the  edge  coefficients  in  defined  in  section  2.1,  ij=l2^  denote  the 
x,y,z  subscripts,  (as  opposed  to  vertex  addresses  as  in  die  2D  case)  and  the  proportionality 
constant  has  been  absorbed  into  the  Oij  coefficients.  The  summation  is  over  all  tetrahedral  ele¬ 
ments  surrounding  the  edge  (i.e  edge  1-2  in  Figure  2).  F,  denotes  the  ith  component  of  the 
normal  area  vector  of  the  face  of  tetrahedron  e ,  which  touches  vertex  2,  but  does  not  contain 
edge  1-2,  and  Gj  the  jth  component  of  the  face  nomial  of  tetrahedron  e  which  touches  vertex 

1,  but  does  not  contain  edge  1-2.  Since  these  coefficients  are  symmetric  in  F  and  G,  they  can 
be  used  to  compute  the  flux  contribution  to  vertex  1  as  well  as  the  flux  contribution  to  vertex 

2,  along  edge  1-2.  This  requires  the  storage  of  9  coefficients  per  edge.  If  the  coefficients  are 
symmetric  in  the  i,j  indices,  then  only  6  coefficients  per  edge  are  required.  In  order  to  prove 
this  property,  we  form  the  difference 


Okj  -  dji 


(12) 


The  volume  of  a  tetrahedral  element  can  be  expressed  as 

vol,  =  (Fi  Gj  -  Fj  Gi)Axt  (13) 

where  i=j=k  and  Ax*  represents  the  difference  (xt)a  -  (x*)*,  where  a  and  b  are  the  two  rwnain- 
ing  vertices  other  flian  the  edge  endpoints  1,2  of  the  tetrahedron  e,  as  shown  in  Figure  2. 
Inserting  this  expression  into  equation  (12),  we  obtain 


z 


Fj  Gj 

(Fi  Gj  -  Fj  G,  )Ax* 


-  I 


Fj  Gj 


^  (  F,  Gj  -  Fj  Gi  )Ax* 


(14) 


Since  the  summation  is  over  all  tetrahedra  surrounding  die  edge  1-2,  the  sum  of  all  such  Axt 
forms  a  closed  polygon  encircling  the  edge,  and  thus  vanishes.  The  coefficients  are  thus  sym¬ 
metric  in  ij  for  all  values  of  i  and  j,  and  thus  only  6  coefficients  per  edge  need  to  be  stored. 
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Flgare  2 

Illustration  of  Tetrahedra  Considered  for  Formuttion  of  Coefficient  Associated  with  Edge  1-2. 
F,  and  G,  rqiresent  Face  Area  Normals 
to  the  Two  Exposed  Faces  cf  Each  Tetrahedron  e . 
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Flgarc  3a 

Fine  Grid  Emirioyed  for  the  Rae  2822  Wing  Test  Case 
(1.04  million  (Kunts,  6  million  leinhedia.  Wall  q»cing  ICT^  chords) 


Figure  3b 

Second  Finest  Grid  Employed  for  the  Rae  2822  Wing  Test  Case 
(135,000  points) 


Figare  3c 

Third  Grid  Employed  for  the  Rae  2822  Wing  Test 
(23,000  pcHnts) 


Figure  3d 

Fourth  Grid  Employed  for  die  Rae  2822  Wing  Test  Case 
(2610  points) 


Figare  4 

Computed  Density  Contours  for  the  Rae  2822  Wing 
(Mach  =  0.73,  Re  =  6.5  million.  Incidence  =  2.79  degrees) 
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Figorc  6 

Comparison  of  Observed  Convergence  Rates  for  the  Three-Dimensional 
Rae  2822  Wing  Multigrid  Case  and  die  Corresponding  Two-Dimensional 
Case  both  widi  and  without  Multigrid  Acceloalion 
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Figare  7 

Computed  Density  Contours  for  die  Swqit  and  Ti^iered  Rae  2822  Wing 
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Figure  S 

Fine  Grid  Employed  for  the  3-Element  Wing  Test  Case 
(1.84  million  points,  10.6  million  tettahedra.  Wall  spacing  l(r®  chords) 


Figare  9 

Ccmpuled  Solution  for  the  3-Eleinent  Wing  Test  Case. 

Mach  contours  are  displayed  on  the  Symmetry  Plane,  and  Density  Contours  on  the  Wing  Surfaces. 
(Mach  =  0.2,  Re  =  9  million.  Incidence  =  1621  degrees) 
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Figure  11 

Computed  Solution  for  the  3-Element  Swept  and  Tapered  Wing  Test  Case. 

Mach  contours  are  displayed  on  the  Symmetry  Plane,  and  Density  Contours  on  the  Wing  Surfaces. 
(Mach  =  0.2,  Re  =  9  n^on.  Incidence  =  12  degrees) 
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